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ABSTRACT 

During their evolution, star clusters undergo mass segregation, by which the orbits of 
the most massive stars shrink, while the lighter stars move outwards from the cluster 
centre. In this context, recent observations and dynamical modelling of several galactic 
and extra-galactic globular clusters (GCs) suggest that most of them show, close to 
their centre, an overabundance of mass whose nature is still matter of debate. For 
instance, many works show that orbitally segregated stars may collide with each other 
in a runaway fashion, leading to the formation of a very massive star or an intermediate 
mass black hole (IMBH) with a mass comparable to the observed mass excess. On the 
other hand, segregated stars can form a dense system if the IMBH formation fails. 
In this paper we study the early formation phase of a dense, massive sub-system 
(MSS) in several GCs models using a recently developed semi-analytical treatment of 
the mass segregation process. In order to investigate how the MSS properties depend 
on the host cluster properties, we varied initial mass function (IMF), total mass, 
spatial distribution and metallicity of our models. Our results show how the IMF 
contributes to determine the final mass of the MSS, while the metallicity and the 
spatial distribution play a minor role. The method presented in this paper allowed 
us to provide scaling relations that connect the MSS mass and the host cluster mass 
in agreement with the observed correlation. In order to follow the early formation 
stage of the MSSs and improve our statistical results, we performed several TV-body 
simulations of stellar clusters with masses between 10 3 and 2 x 10 5 solar masses. 

Key words: stars: stellar evolution; stars: black holes; stars: kinematics and dynam¬ 
ics; galaxies: star clusters; Galaxy: globular cluster. 


1 INTRODUCTION 

Recently, a number of observations of globular clusters 
(GCs) suggested that their cores may contain much more 
mass then expected from earlier observations (van der Marel 
& Anderson 2010; Noyola et al. 2010). 

As proposed by several authors, these mass excesses can 
be ascribed to the presence of an intermediate mass black 
hole (IMBH) or a very massive star (VMS) close to the clus¬ 
ter centre, with masses in the range 10 2 — 10 4 Mg . 

For instance, observations and modelling of the inner¬ 
most region of M15 seem to be consistent with the pres¬ 
ence of a single object (Gerssen et al. 2002; den Brok et al. 
2014). Furthermore, hints for IMBHs candidates have been 
found also in other GCs in the Milky Way (Noyola et al. 
2008; Liitzgendorf et al. 2013; Feldmeier et al. 2013; Peuten 
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et al. 2014), and in other galaxies (see for example Mapelli 
et al. (2008)), as in the case of G1 in M31 (Gebhardt et al. 
2002, 2005; Miller-Jones et al. 2012), and the young mas¬ 
sive cluster in the M82 irregular galaxy (Kaaret et al. 2001; 
Matsumoto et al. 2001; Usuda et al. 2001). 

Several mechanisms have been proposed for the birth 
and growth of IMBHs. For instance, Miller & Hamilton 
(2002) pointed out that a stellar BH seed may grow slowly 
through occasional collisions and merging with other stars. 

Another possibility relies upon the mass segregation 
process, driven mainly by dynamical friction (df), which 
leads to an accumulation of mass within the cluster cen¬ 
tre in form of orbitally decayed stars. In such a case, the 
most massive stars tend to concentrate toward the centre of 
the cluster, while the lighter component moves outward in 
an attempt to establish energy equipartition. As the mass 
segregation proceeds, the heaviest stars lose kinetic energy 
and reach the innermost region of the cluster, where they 
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form a dense, contracting nucleus (Spitzer 1969; Heggie & 
Hut 2003; Binney & Tremaine 2008). 

In this framework, some authors suggested that the 
shrinking nucleus collapses in a fraction of the relaxation 
time-scale (Portegies Zwart & McMillan 2002; Giirkan et al. 
2004). The collapse facilitates a phase of runaway collisions 
among stars that leads to the formation of an IMBH (Porte¬ 
gies Zwart et al. 1999, 2004; Freitag et al. 2006; Goswami 
et al. 2012; Liitzgendorf et al. 2015; Giersz et al. 2015), 
whose gravitational field can significantly shape the dynam¬ 
ical evolution of the host cluster (Liitzgendorf et al. 2013; 
Leigh et al. 2014). 

Another possibility is that the observational evidence 
of a compact object in the centre of GCs can be interpreted 
as a dense, massive sub-system composed of dark remnants 
of heavy stars (Baumgardt et al. 2003; van der Marel & 
Anderson 2010; Haggard et al. 2013; Lanzoni et al. 2013; 
Kamann et al. 2014). In this case, the formation of binaries 
efficiently halts the contraction of the nucleus, preventing 
the formation of a VMS (or an IMBH). The evolution of 
the resulting sub-system of stars, which we refer to as MSS, 
will be likely dominated by two and three-body interactions 
(Baumgardt et al. 2003; Trani et al. 2014). 

A great effort toward the understanding of IMBHs for¬ 
mation and growth will likely come from the next generation 
of space-based gravitational waves (GWs) observatories, as 
the eLISA satellite (Amaro-Seoane & Freitag 2006; Amaro- 
Seoane et al. 2009; Mapelli et al. 2010) but, currently, it is 
quite hard to discriminate between an IMBH and an MSS 
through observations. 

In this paper, we investigate how the global properties 
of the host cluster affect the formation process of an MSS. 
Using a treatment for the dynamical friction process recently 
developed by Arca-Sedda & Capuzzo-Dolcetta (2014a), we 
studied mass segregation in 168 models of GCs with differ¬ 
ent total masses, density distributions, initial mass functions 
(IMFs) and metallicities. 

Our results suggest that the kind of stars which pop¬ 
ulate the MSS depends on the IMF and metallicity of the 
host cluster, while spatial distribution and stellar evolution 
affect significantly the MSS mass. Furthermore, we provide 
relations connecting the MSS mass with the total mass of 
the host cluster, showing that the best agreement with ob¬ 
servations is achieved for power-law IMFs. 

To follow the MSS formation process in detail, we 
performed direct A-body simulations of several GCs with 
masses in the range 10 3 — 10 s Mg. The time-scale needed to 
assembly an MSS, as well as the MSS mass and size, agrees 
with our semi-analytical results. 

In Section 2 the methodology used to derive MSS 
masses is introduced and discussed, in Section 3 we inves¬ 
tigate the properties of MSSs and provide scaling relations 
connecting the mass of MSSs and the host cluster, in Sec¬ 
tion 4 we discuss the direct IV-body simulations performed 
whereas in Section 5 we draw the conclusions of this work. 


2 NUMERICAL METHOD 


with them (see for example van der Marel & Anderson 
(2010); Haggard et al. (2013); Lanzoni et al. (2013)). We 
will discuss the implications of different observational mass 
estimates on our results in the next section. Table 1 sum¬ 
marises the main properties of the observed GCs. 

To reach our aims, we need two important ingredients: 
i) a detailed description of the df process, which primarily 
drives heavy stars toward the host cluster centre, and ii) a 
reliable modelling of the host GC. 

In order to estimate the final mass of MSSs, we sampled 
several GC models by varying their total mass, radius, initial 
mass function (IMF) and metallicity (Z). In particular, for 
each model we selected initial position, orbital eccentricity 
and mass of all its stars, which are crucial ingredients to 
determine how much stars contribute to the formation of an 
MSS. 


2.1 Dynamical friction 

A massive body traversing a sea of lighter particles suffers 
a dynamical braking that drags it toward the centre of the 
host system (Chandrasekhar & von Neumann 1943; Chan¬ 
drasekhar 1943a,b). 

This mechanism, called dynamical friction (df), arises 
directly from two-body encounters and plays a signifi¬ 
cant role in shaping the evolution of astrophysical sys¬ 
tems on very different scales (Bekenstein 1973; Tremaine 
1976; Capuzzo-Dolcetta 1993; Milosavljevic & Merritt 2001; 
Gualandris & Merritt 2008; Antonini 2013; Arca-Sedda & 
Capuzzo-Dolcetta 2014b; Arca-Sedda et al. 2015). 

In a pioneriing paper, Chandrasekhar (1943a) provided 
the timescale over which a body of mass m* that moves 
on an orbit with initial apocenter r* and initial velocity n», 
reaches the innermost region of its host system (Binney & 
Tremaine 2008): 


rdf(Myr) 


1.9 x 10 4 / r, \ 2 / v, \ / 10 8 M 0 \ 

logA \5 kpc ) \200 kms -1 / \ m* ) ’ 

( 1 ) 


where, for star clusters, we assume logA ~ 10 for the usual 
Coulomb logarithm. 

Many works attempted to generalise the Chan¬ 
drasekhar’s work to axisymmetric and triaxial systems (Bin¬ 
ney 1977; Ostriker et al. 1989; Pesce et al. 1992), and to sys¬ 
tems characterised by cusped density profiles (Merritt et al. 
2006; Vicari et al. 2007; Antonini & Merritt 2012; Arca- 
Sedda & Capuzzo-Dolcetta 2014a). 

In particular, Arca-Sedda & Capuzzo-Dolcetta (2014a) 
developed a reliable treatment of df particularly well suited 
to describe the motion of massive bodies in both cusped and 
cored density profiles. Furthermore, they provided an use¬ 
ful formula for the df time-scale, recently updated in Arca- 
Sedda et al. (2015): 


rdf(Myr) = T re f 


r 3 
' GC 

Mac 


g(e, 7 ) 


Mac J 


\ - 0.67 / 

\r gc 


_ ( 2 ) 

where T re f = 0.3\/ 10 11 MQ/lkpc 3 and e identifies the 
initial eccentricity of the orbit 


As database to compare with our theoretical results, we used 
observational data provided by Liitzgendorf et al. (2013) 
(hereafter LU13), although some works seem to be at odds 
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with r v the pericentral distance from the cluster centre. The 
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Table 1. Parameters of the observed GCs collected in LU13. 


ID 

NAME 

Log(M GC /M 0 ) 

£ M tot 

Log(M B H/M Q ) 

e -M BH 

G1 


6.76 

0.02 

4.25 

0.11 

NGC104 

47Tuc 

6.04 

0.02 

< 3.17 

— 

NGC1851 


5.57 

0.04 

< 3.3 

- 

NGC1904 

M79 

5.15 

0.03 

3.47 

0.12 

NGC2808 


5.91 

0.04 

< 4 

- 

NGC5139 

cjCen 

6.40 

0.05 

4.6 

0.08 

NGC5286 


5.45 

0.02 

3.17 

0.24 

NGC5694 


5.41 

0.05 

< 3.9 

- 

NGC5824 


5.65 

0.03 

< 3.78 

— 

NGC6093 

Mso 

5.53 

0.03 

< 2.9 

- 

NGC6266 

M62 

5.97 

0.01 

3.3 

0.18 

NGC6388 


6.04 

0.08 

4.23 

0.18 

NGC6715 

M54 

6.28 

0.05 

3.97 

0.18 

NGC7078 

M15 

5.79 

0.02 

< 3.64 

- 


Column 1: name of the cluster. Column 2: alternative name of the cluster. Column 3: mass of the cluster. Column 4: logarithmic error 
on the cluster mass. Column 4: mass of the central IMBH candidate. Column 5: logarithmic error on the IMBH candidate mass. 


function g(e), instead, links e and the slope of the density 
profile, 7 , through the relation 


fl(e, 7) 


(2-7) 


ai 


((2-7 )“ 2 



(1 


e) +e 


(4) 


with a\ = 2.63 ± 0.17, a 2 = 2.26 ± 0.08 and 03 = 0.9 ± 0.1. 

Considering a typical GC characterised by a cored den¬ 
sity profile (7 = 0 ), with a mass Mgc ~ 10 6 Mq and a scale 
radius roc = 1 pc, it easy to find through Equation 2 


Tdf < 0.6 Gyr, (5) 

for a massive star (m* ~ 25 Mq) that moves on a circular 
orbit at r* = 5 pc from the GC centre. 

Hence, a population of heavy stars may sink to the cen¬ 
tre of the host GC in few Gyr, leading to a significant accu¬ 
mulation of mass in its innermost region. 


2.2 Sampling method 

To provide reliable estimates of the amount of stars that 
have sunk to the cluster centre in a given time, we selected 
isolated cluster models with masses in the range 10 3 — 3 x 
10 6 Mq, which differ each one in spatial distribution, IMF 
and metallicity. Furthermore, we included the SSE package 
(Hurley et al. 2000) in our statistical code, in order to take 
in account stellar evolution, which causes stellar mass loss 
and may alter the df time. The properties of our GC models 
are discussed in detail in the following sections. 


Spatial distribution 

We sampled the positions of each star according either to an 
uniform spatial distribution, in which the density of stars 
does not depend on the spatial coordinates, or to a cored 
7 -profile (7 = 0) (Dehnen 1993). 

Furthermore, we assigned to stars’ orbits an orbital ec¬ 
centricity, e, in the range e = 0 (circular orbits) and e = 1 
(pure radial orbits) according to a flat distribution. We ver¬ 
ified that the use of different sampling methods for e does 
not affect significantly the global results, unless all the stars 
move on circular, or purely radial, orbits. 


Although an uniform sphere evolves over few dynami¬ 
cal times and the treatment presented here do not account 
for the time evolution of its global structure, we will show 
in the next section that our semi-analytical estimates well 
reproduce the accumulation of mass in cluster models char¬ 
acterised by a constant density by comparing them with 
direct IV-body simulations. 


Initial mass function and metallicity of the cluster 

We assigned masses to the stars in the range 0.1 and 100 
Mq according to either a flat IMF, a Salpeter IMF (Salpeter 
1955) or a Kroupa IMF (Kroupa 2001). 

Moreover, to highlight the effects of different metallici- 
ties (Z) on the formation of an MSS, we assigned to stars 
either a solar metallicity Z = 0.02, or the typical metallicity 
of old globular clusters, Z = 0.0004. 

At the end, we gathered a total sample of 168 models, 
whose main parameters are summarised in Table 2. 

Moreover, we made 100 realisations of each cluster, in 
order to filter out statistic fluctuations, providing an aver¬ 
aged value of the MSS mass along with an estimate of the 
corresponding standard deviation. 

We grouped clusters having different masses but the 
same global properties, labelling them with the letter A or 
B and a number between 1 and 6 . 

It should be noted that such a choice brings together 
systems in which the typical time-scale over which an MSS 
form can be very different, since it is directly connected to 
the total mass of the host cluster. 

In particular, our grouping choice gathers clusters with 
ages greater than the typical time-scale over which an MSS 
forms. 

Letter A refers to models with solar metallicity, whereas 
letter B indicates metal-poor models. The number, instead, 
depends on the kind of IMF and the spatial distribution of 
the stars. Each number represents a combination of the IMF 
and the spatial distribution that characterises the model. For 
instance, number 1 refers to models with a Kroupa IMF and 
an uniform spatial distribution. 
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Table 2. 


Parameters of the cluster models. 


Model 

IMF 

p(r) 

Z 

Mgc ( m q) 

Al 

K 

U 

0.02 

10 3 

- 3 x 10° 

A2 

K 

D 

0.02 

10 3 

- 3 x 10 6 

A3 

S 

U 

0.02 

10 3 

- 3 x 10 6 

A4 

S 

D 

0.02 

10 3 

- 3 x 10 6 

A5 

F 

U 

0.02 

10 3 

- 3 x 10 6 

A6 

F 

D 

0.02 

10 3 

- 3 x 10 6 

B1 

K 

U 

0.0004 

10 3 

- 3 x 10 6 

B2 

K 

D 

0.0004 

10 3 

- 3 x 10 6 

B3 

S 

U 

0.0004 

10 3 

- 3 x 10 6 

B4 

S 

D 

0.0004 

10 3 

- 3 x 10 6 

B5 

F 

U 

0.0004 

10 3 

- 3 X 10 6 

B6 

F 

D 

0.0004 

10 3 

- 3 X 10 6 


Note. Column 1: name of the model. Column 2: IMF used to 
sample masses of the stars: Kroupa (K), Salpeter (S) and flat 
(F). Column 3: spatial density profile of the stars: uniform (U) 
and Dehnen (D). Column 4: metallicity of the cluster. Column 5: 
masses of the cluster models. 

As pointed out above, we selected three IMFs: flat (de¬ 
noted with letter F in Table 2), Kroupa (K) and Salpeter 
(S). The spatial distributions considered, instead, are uni¬ 
form (U) or cored 7 -distribution (D). 


2.3 Evaluation of MSS masses and sizes 

A star undergoes orbital decay only if its mass exceeds the 
mean mass of the background stars (m*) (Chandrasekhar 
1943a). In particular, several works pointed out that dy¬ 
namical friction affects significantly the motion of bodies 
with masses m* ^ K(m *), where A is a parameter in the 
range 1 — 50, depending on the method used in describing df 
(Colpi et al. 1999; Arca-Sedda & Capuzzo-Dolcetta 2014a; 
Alessandrini et al. 2014; Miocchi et al. 2015). In this paper, 
we assumed K = 30, a value well supported by previous 
theoretical and numerical results (Arca-Sedda & Capuzzo- 
Dolcetta 2014a). 

Furthermore, mass loss process can alter the decay pro¬ 
cess, leading to a significant increase of the df time, at least 
for the most massive stars. 

In particular, it is worth nothing that stars with masses 
above 30 Mg lose most of their initial mass in < 10 Myr, 
leading to a significant decrease of the df efficiency well be¬ 
fore they reach the cluster centre. Indeed, an inversion of 
Equation 2 allows to guess the initial position that a star 
should have to reach the cluster centre within a time t: 

r(t,m») = (^=) (£,) roc, (6) 

where T = T re fp(e, 7 ) (r GG /M GG ) 1/2 ■ 

For example, a star with m, = 30 Mg that moves in a 
typical GC with M GG = 10 6 Mg, r GG = 1 pc and 7 = 0, 
should have an initial position r(t , m*) ~ 0.6 pc to reach the 
GC centre within 10 Myr. 

Hence, we took in account how the df time-scale changes 
as a consequence of the stellar evolution proceeding in the 
following way. We divided the time t in sub-intervals 5ti 


equally spaced such that JA 5ti = t, keeping each interval 
small enough to follow in detail the evolution of the star mass 
as a function of the time, i.e. 5ti < 1 Myr. For each value of 
5ti, we evaluated the updated mass of the star through the 
SSE package, m*;, and the df time that corresponds to such 
a mass, Tdf,i, using Equation 2. 

We defined the final df time-scale of a star as the 
weighted-average value of all the df times 5ti over the to¬ 
tal time interval: 


bi r 


^ ~2j m,jTdf,i 
Ei m *i 


( 7 ) 


using its actual mass m*i as weight. 

At the end, the total mass of the MSS is then evaluated 
as the sum of the masses of those stars with tdf C t and 
m* > K{m), plus the masses of stars with initial position 
smaller than the expected MSS size, r* < rMss- 

Furthermore, this method allows to provide hints on the 
size of the MSS. Indeed, as pointed out by several authors 
(Kalnajs 1972; Read et al. 2006; Gualandris & Merritt 2008; 
Antonini & Merritt 2012; Arca-Sedda & Capuzzo-Dolcetta 
2014a), df stops when the star, with mass m*, approaches a 
distance from the cluster centre that encloses an amount of 
mass roughly equal to m*. This critical radius is commonly 
called “stalling radius”, rMss- Within rMSS, the motion of 
the star is mainly dominated by random encounters, while 
the df action becomes negligible. Since the greater the star 
mass the greater the stalling radius, we can define as MSS 
size the tmss of the most massive star (m* = 100 Mg). 

For a 7 -profile, whose radial mass distribution is given 
by: 

/ \ 3 — 7 

M GG (r) = Mq G ( — A- - , ( 8 ) 

\r + r gc / 


this “stalling radius” can be obtained by inverting the latter 
equation, leading to: 

M , Q , 

mss=7'Gc^—(9) 

with M = (t-tu/Mgc) 1 ^ 3 ^■ 

For an uniform distribution, instead, it is easy to find: 

r M ss = Rgc ( 4 ^) , (10) 


being R G c the total radius of the cluster. 

Figure 1 shows the MSS sizes as a function of the GC 
mass for different values of the cluster scale radius, r GG , in 
the case of a 7 density profile and in the case of an uniform 
density profile. 

It is quite evident that the heavier the host cluster the 
smaller the MSS radius, which can reach values below 0.1 
pc for Mq G > 10 6 Mg. It is worth noting that most of the 
observed mass excesses are enclosed within a region whose 
extension is < 0.5 pc, quite close to the radii of our MSSs 
(van der Marel & Anderson 2010; Haggard et al. 2013; Lan- 
zoni et al. 2013). 

The time-scale over which the formation of an MSS 
takes place depends on the total mass of the host cluster. In 
particular, as the stars’ orbits shrink to the cluster centre, 
the MSS mass increases, reaching a saturation value over a 
time-scale that represents an upper limit to the relaxation 
time of the system (Binney & Tremaine 2008; Henon 1960; 
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Figure 1. Size of the MSSs obtained using Equation 9 for models 
with a 7 density profile (top panel) and Equation 10 for models 
with an uniform density profile (bottom panel). 

Freitag et al. 2006). Hence, it is evident that our choice to 
group models with the same properties but different masses 
implies that the time-interval over which their MSS should 
have formed can be very different. 

Though our methodology is substantially blind to their 
subsequent evolution, we examine in the following two pos¬ 
sible fates for the MSSs for the sake of completeness: 

• heavy stars accumulate into the cluster centre in a time 
much shorter than the stellar mass loss time-scale. In this 
case, the contraction of the MSS drives the inner region of 
the cluster toward core collapse, thus facilitating a runaway 
collision phase (Portegies Zwart et al. 1999; Portegies Zwart 
& McMillan 2002; Portegies Zwart et al. 2004; Freitag et al. 
2006; Liitzgendorf et al. 2015), even boosted by binaries 
formed during the collapse process (Portegies Zwart et al. 
1999; Portegies Zwart & McMillan 2002; Freitag et al. 2006; 
Gaburov et al. 2008) and by primordial binaries (Portegies 
Zwart & McMillan 2007), which leads to the formation of 
an IMBH (Portegies Zwart et al. 1999, 2004; Giirkan et al. 
2004). In this scenario, the MSS can represent the total 
reservoir of mass available in form of stars to build-up the 
IMBH through runaway collisions; 

• stellar mass loss occurs over a time comparable to 
the MSS formation. In this case, the formation of binaries 
and multiple systems within the MSS halts efficiently the 
contraction process, thus quenching stellar collisions (An- 
geletti & Giannone 1977; Chernoff & Weinberg 1990; Porte¬ 


gies Zwart & McMillan 2007; Vesperini et al. 2009; Lamers 
et al. 2013; Mapelli & Bressan 2013; Fujii & Portegies Zwart 
2014). In this case, the inner region of the cluster undergoes 
a series of contraction and re-expansions called gravother- 
mal oscillations (Bettwieser & Sugimoto 1984; Cohn et al. 
1989; Makino 1996) that drives an expansion of the MSS by 
a factor up to ~ 2.5 (Trani et al. 2014). It is worth noting 
that our estimates of the MSSs mass represent quite well 
the expected mass of these small cores. 


3 RESULTS 

Using the approach described above, we examine in this sec¬ 
tion the stellar content of the newly born MSSs, providing 
also correlations that link the masses of MSSs and those of 
their host clusters. 

3.1 Properties of MSSs 

In the following, we use models A2 and B2 as reference cases, 
in order to highlight the differences arising from the choice of 
different metallicities, and models Al and A2 to highlight, 
instead, the effects connected to different spatial distribu¬ 
tions of the host cluster. As pointed out in the previous 
sections, we would demonstrate here that the mass excess 
observed in the centre of several GCs is likely related to a 
sub-system of orbitally segregated stars, not necessarily to 
an IMBH. 

Figure 2 shows the number of orbitally segregated stars, 
Nmss, as a function of the time for models Al, A2 and B2 
and for different GC total masses. 

The growth process is characterised by two distinct 
phases: one more rapid, lasting up to 0.1 — 0.5 Gyr, during 
which stars with nearly radial orbits or with small apoc- 
entres segregate fastly to the centre, and a second, slower 
phase, during which the main contribution to the MSS 
growth is given by stars that move on a more peripheral 
region of the cluster, thus reaching occasionally the GC cen¬ 
tre. 

However, since stars lose mass during their evolution, 
the parameter Nmss cannot be used to provide informations 
about the MSS properties. In order to investigate how the 
MSS mass increases in time, we show in Figure 3 the total 
amount of mass accumulated within the cluster centre for 
the same models, considering a GC model with mass M = 
10 6 M 0 . 

Looking at the figure, three different phases are clear: 
the deposited mass increases until it reaches a maximum 
value, then the mass accumulation process undergoes a rapid 
decrease phase which last 2 — 3 Myr and finally it rises 
smoothly toward a saturation value. The first phase cor¬ 
responds to the decay of the most massive stars located in 
an inner region of the GC. Then, the mass decrease during 
the second phase as a consequence of the mass lost by segre¬ 
gated stars. Finally, the last stage is mostly determined by 
the deposit of stars with initial orbits located in the outer 
shells of the GC. 

Comparing the MSS mass growth for model Al and 
A2 in the case of Mgc = 10 6 M@, it is evident that the 
deposited mass in model A2 is initially greater than in model 
Al. This is due to the fact that the A2 density profile is more 
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Log t (Myr) 

Figure 2. Number of orbitally segregated stars as a function of 
the time for different masses of the host GC. Each panel refers 
to a different model, as indicated. Within each panel, instead, 
from bottom to top, lines refer to a cluster model with a mass 
M(jc = 10 4 — 10* — 10® — 3 X 10® Mq, respectively. 


concentrated, thus implying that in this case stars move on 
orbits with smaller apocentres (on average) and therefore 
smaller df times. On the other hand, over a time > 2 Gyr, 
the deposited mass reach comparable values both in model 
A1 and A2. 

This is mainly due to the fact that the density profile 
of model A2 scales as p(r) oc (r + tgc) 4 and, therefore, 
stars with initial apocenter r r gc travel in a low dense 
environment in which df is highly suppressed. On the other 
hand, since model A1 has an uniform density profile, stars 



Log t (Myr) 



Figure 3. Comparison between the mass deposited in the centre 
of a cluster with M = 10® Mq for models A1 and A2 (top panel) 
and for models B2 and A2 (bottom panel), respectively. 

moving on farther orbits may reach the GC centre since df 
acts efficiently also on larger length-scales. 

Comparing models A2 and B2, instead, we found that 
the MSS mass evolution is very similar. In particular, in 
model B2 the mass of the MSS reaches a saturation value 
greater than in model A2. This is due to the fact that stars 
with masses above 10 Mq having low metallicities reach final 
masses up to 15% greater than stars with solar metallicity. 
Such a difference may be even greater depending on the 
kind of stellar evolution recipes considered, as pointed out 
by several authors (Brocato et al. 1999; Mapelli & Bressan 
2013; Ziosi et al. 2014; Spera et al. 2015). 

We can also provide hints about the stellar content of 
the MSSs, looking at stars form them. In particular, Figure 
4 shows the fractional number of stars that should form the 
MSS of a GC with a total mass Mgc = 10® Mq and an 
age t ~ 6 Gyr. Stellar types are defined as in Hurley et al. 
(2000), and are listed in Table 3. 

The comparison between models with different metal¬ 
licities (A2 and B2) highlights that a cluster with lower val¬ 
ues of Z will host an MSS likely dominated by BHs, which 
are more than 60% of the total number of orbitally segre¬ 
gated stars, while neutron stars (NSs) seem to be the most 
common stars in a MSS of a cluster with solar metallicity. 
Furthermore, also the spatial distribution of stars seems to 
be important in determining the stellar composition of the 
MSS. Indeed, the fraction of BHs and NSs is smaller in mod¬ 
els with an uniform spatial distribution (Al) with respect to 
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Table 3. 


Evolution phases. 


REF N. 

NAME 

Description 

0 

MSI 

MS stars with M ^ 0.7Mq 

1 

MSh 

MS stars with M > 0.7M© 

2 

HG 

Hertsprung Gap 

3 

GB 

First GB 

4 

CHeB 

Core Helium Burning 

5 

EAGB 

Early AGB 

6 

TPAGB 

Thermally Pulsing AGB 

7 

HcMS 

Naked Helium star MS 

8 

HcHG 

Naked Helium star Hertsprung Gap 

9 

HeGB 

Naked Helium star GB 

10 

HcWD 

Helium White Dwarf 

11 

COWD 

Carbon/Oxygen White Dwarf 

12 

ONeWD 

Oxygen/Neon White Dwarf 

13 

NS 

Neutron Star 

14 

BH 

Black Hole 

15 

mO 

massless remnant 


40 

30 

20 


10 



Figure 4. Number of orbit ally decayed stars, within t = 6 Gyr, 
of a given stellar type over the total number of decayed stars in 
models A2 and A1 (top panel) and in models A2 and B2 (bottom 
panel) for a cluster with M = 10 6 Mq. 


the models with a 7 density profile (A2). Regarding model 
Al, we found that also a population of low main sequence 
stars (MSI) with inital apocentres quite close to the centre 
of the host GC, can contribute significantly to the total MSS 
mass. 

Equation 2 can be used to provide some hints about the 
global time evolution of the stellar distribution of the stars 
in the cluster. In particular, the position, r, of a star at a 
given time t can be obtained by solving the equation: 

Tdf(r*) - rdf(r) = t, (11) 


whose explicit solution is given by: 


r(t) = r . 



Tdf(r*) 


( 12 ) 


Hence, Equation 12 can be used to investigate how the 
spatial distribution of stars evolves in time. 

Figure 5 shows the time evolution of the half-mass ra¬ 
dius rh of main sequence stars (MSs), white dwarfs (WDs), 
neutron stars (NSs) and black holes (BHs) for a GC with 
M = 10 6 Mq in models Al, A2 and B2. The interval of time 
considered is comparable to the time-scale needed to build 
up an MSS, as shown in Figure 3. It is evident that in all 
the cases considered, the population of BHs and NSs tends 
to concentrate in an inner region of the host GC, while WDs 
and MSs stars move in an outer region, having their motion 
less affected by df. 


Column 1: reference number used in Hurley et al. (2000). Column 
2: abbreviation. Column 3: stellar evolution phase. 


As shown in Figure 6 , the population of BHs concen¬ 
trate more efficiently than other stellar types, reaching val¬ 
ues rh,BH < 0.1 pc quite independently from the spatial 
distribution or the metallicity of the host cluster. 

Our results suggest that MSSs are mainly composed of 
BHs and NSs. On the other hand, it is currently quite ascer¬ 
tained that NSs can receive a kick at their birth caused by 
possible asymmetries that occur during the core-collapse of 
the NS progenitor (Gott et al. 1970; Harrison & Tademaru 
1975; Lai et al. 2001). Though the dispersion and distribu¬ 
tion of the kick velocities is quite uncertain, most of the kicks 
estimates derived on theoretical and observational basis can 
throw the NS away from its host cluster (see for example 
Podsiadlowski et al. (2005) for a general discussion on this 
topic). The loss of NSs in the cluster can have significant 
effects on its global long term evolution, as shown recently 
by Contenta et al. (2015). 

Kick velocities are widely described through maxwellian 
distributions with a in the range 190—260 km s -1 (Hansen & 
Phinney 1997; Hobbs et al. 2005). However, this description 
seems to be at odds with observations, at least in part. In¬ 
deed, several GCs in the Milky Way contain more than 10 3 
NSs, whereas the kick velocities cited above should leave 
in the host cluster few tens of them. This problem, com¬ 
monly referred to as “NS retention problem”, can be par¬ 
tially solved if a substantial fraction of NSs form in massive 
binaries, where the kick momentum suffered by the NSs at 
their birth is shared with the companion, lowering the net 
velocity of the system (Drukier 1996; Ivanova et al. 2005; Da- 
lessandro et al. 2011). Some authors proposed that some NSs 
receive low birth kicks (< 50 km s -1 ) depending on the prop¬ 
erties of the binary in which they form (Pfahl et al. 2002; 
Podsiadlowski et al. 2004), or depending on the amount of 
mass that fall back during the core collapse of the star (Fryer 
et al. 2012 ). 

In order to take in account which effect natal kicks can 
have on the population of NSs that compose the MSSs, we 
estimate in the following the fraction of NSs that are ejected 
in our GC models assuming that the velocity kicks are de- 
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Figure 5. Half-mass radius of different types of stars as a function 
of the time for a cluster with total mass Mqc = 10® Mq. From 
top to bottom, panels refer to model Al, A2 and B2, respectively. 

scribed by a maxwellian with a = 190 km s -1 as suggested 
by Hansen & Phinney (1997). In particular, we restrict this 
study applying natal kicks to models characterised by a 
Kroupa distribution of stellar masses and a Dehnen density 
profile. 

It is trivial to show that the local escape velocity from 
a given position ro is defined as 

iw(r 0 ) = (20(r o )) 1/2 , (13) 

where <p(ro) represents the gravitational potential generated 
by the cluster evaluated at the position ro. Once the NS 
is kicked out, dynamical friction erases part of the kinetic 



t (Gyr) 



t (Gyr) 

Figure 6. Half-mass radius of BHs as a function of the time 
for a cluster with total mass Mq c = 10® Mq. Top panel refers 
to models Al (straight line) and A2 (dotted line) while bottom 
panel refers to models B2 (straight line) and A2 (dotted line), 
respectively. 


energy excess generated by the kick. Indicating the energy 
depleted through df as -Edf(fo), the correct escape velocity 
is given by: 

Vcor(ro) = fesc(ro) ^1 + ~^ry ^ • ( 14 ) 

Since 0(0) = {GMgc/tgc) for a Dehnen model, the 
escape velocity from the innermost region of a GC with 
Mgq = 10 6 Mq and tgc = 0.15 pc, is it?sc(0) ~ 170 km 
s -1 . 

Following the numerical approach described by Arca- 
Sedda & Capuzzo-Dolcetta (2014a), we found that the term 
■E'dt(fo)/</ , (*"o) never exceeds 10 -2 for a typical kick velocity 
in the range 10 2 — 10 3 km s _1 , thus implying an increase of 
the escape velocity at most of 1%. 

Therefore, assuming a maxwellian distribution of the 
kicks with a = 190 km s~ 3 and an escape velocity of ~ 170 
km s -1 , we found that ~ 10% of the NSs that compose the 
MSS should remain bound to the host cluster, and therefore 
can come back to the MSS over a df time-scale. However, this 
fraction can be much lower, depending on the GC properties 
and the dispersion of the velocity kicks, as shown in Table 
4. 
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Table 4. 

Escape velocity for a Mqq = 10 6 Mg. 


r GC 

pc 

^cor(0) 

km s —1 

/260 

% 

/l90 

% 

0.1 

230 

15 

31 

0.2 

148 

4.5 

10.5 

0.5 

93 

1.2 

2.9 

1.0 

66 

0.4 

1 


Column 1: scale radius of the GC. Column 2: escape velocity 
evaluated through Equation 14. Column 3: percentage of retained 
NSs assuming a = 260 kms -1 . Column 4: percentage of retained 
NSs assuming a = 190 kms -1 . 

3.2 Scaling relations 

In this section we provide scaling relations that link the MSS 
mass and the host GC mass. 

As database to compare with our results we used the ob¬ 
servational mass estimates for 14 putative IMBHs provided 
by LU13. In particular, LU13 have shown that the mass of 
the IMBH candidates is connected to the host cluster mass 
through the relation 

Log (^) - (101±0.34)Log (Igf) - (2.37±0.34). 

(15) 

However, it should be noted that this relation is charac¬ 
terised by a relative error on its fitting parameters that ex¬ 
ceeds 30%, thus implying quite large errors on the evaluation 
of the IMBH mass for a given mass of the host cluster. Such 
errors are mostly determined by the relatively small amount 
of data available and the large errors that characterise them. 
Hence, here we compare only the statistical behaviour of the 
scaling relation, whereas a more robust comparison need a 
larger amount of data and a much higher precision in the 
evaluation of the mass excess. 

Looking at Figures 7 and 8, which show the MSS masses 
versus the initial mass of the host GCs for all the cases con¬ 
sidered, it is evident that the best agreement with the obser¬ 
vational estimates provided by LU13 is achieved for models 
characterised by a Kroupa IMF, labelled with numbers 1 
and 2, and a Salpeter IMF, labelled with numbers 3 and 4. 

As expected, models with a flat IMF, labelled with num¬ 
ber 5 and 6, produce MSSs with masses that are clearly much 
smaller than the observational estimates. 

Nevertheless, a flat IMF is completely at odds with ob¬ 
servations, and therefore it is not surprising to find MSS 
masses very different from observational estimates in these 
models. 

For all the cases investigated here, we found the follow¬ 
ing power-law best-fitting formula: 

< w) 

where the slope, o, and the offset, b, of the correlation were 
computed by means of a Marquardt-Levenberg non-linear 
regression algorithm. The values of a and b are summarised 
in Table 5. 

It is worth noting that different values of Z produce 
small differences in the correlations. In particular, we found 
that models with low metallicity produces MSSs with masses 
5% greater than that obtained for models with solar metal- 


licities, thus implying at most a difference smaller than 1% 
on the determination of the parameters of the best-fitting 
formula. 

Comparing our best-fitting parameters with that pro¬ 
vided in LU13, it is evident an overall agreement, despite 
the LU13 results are affected by uncertainties that in some 
cases exceed 2 dex. 

Some estimates contained in the LU13 dataset seem to 
be at odds with other works. For instance, some observations 
suggest that the cluster lj Cen hosts in its innermost region 
a mass in the range 1.2 — 1.8 x 10 4 Mg (see for example van 
der Marel & Anderson (2010) and Haggard et al. (2013)), 
less than a half with respect to the value suggested by LU13. 
Another interesting case is the cluster NGC6388, for which 
Lanzoni et al. (2013) suggested a central mass of 2x 10 3 Mg, 
whereas LU13 found a greater value, ~ 1.7 x 10 4 . Using 
the scaling relations drawn here, we can give an estimate 
of the central mass hosted in these two systems. To model 
the old, massive cluster oo Cen, we assume a low value of 
the metallicity, a Kroupa IMF and a cored density profile. 
Note that this choice corresponds to models labelled with 
B2. Using the corresponding parameters (see Table 5) and 
assuming a total mass Mgc = 2.5 x 10 6 Mg (see Table 1), we 
found Mmss = (1.45±0.03) x 10 3 Mg, in agreement with the 
values provided by van der Marel & Anderson (2010) and 
Haggard et al. (2013). Applying the same procedure to the 
cluster NGC6388, we found Mmss — (6.4 ± 1.2) x 10 3 Mg, 
a value which lies between the estimates provided by LU13 
and Lanzoni et al. (2013). A numerical study focused on a 
detailed modelling of this cluster may help to shed light on 
the possible causes of these discrepancies (Liitzgendorf et al. 
2015). 

Finally, for the sake of comparison, we investigated 
whether our results agree with previous theoretical works. 
For instance, using Monte Carlo simulations, Giirkan et al. 
(2004) have shown that the mass of the shrinking core that 
forms in GCs as a consequence of mass segregation, is linked 
to the cluster mass through the relation: 

^1?He(tf)- (2 - s±o - 2) ' (i7) 

in quite good agreement with our estimates both for solar 
and low metallicities. 

Portegies Zwart & McMillan (2002), instead, found that 
the maximum mass of an IMBH formed through runaway 
collisions in a GC with a mass above 10 6 Mg is given by: 

The similarity between our scaling relations and Equa¬ 
tion 18 confirm that MSSs represent the ideal reservoir of 
stars which can contribute to the assembly of an IMBH in 
a runaway fashion. 


4 JV-BODY MODELLING OF A STELLAR 
CLUSTER: IMF, MASS SEGREGATION 
AND MSS EARLY FORMATION 

In order to compare with our statistical results, we ran 10 
direct A-body simulations of star clusters characterised by 
broad IMFs and masses up to 3 x 10 5 Mg. In particular, 
these models consist of 8 — 16 — 32 — 262 — 524k particles. 
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Figure 7. Segregated mass vs. the total mass of the cluster. Black filled circles represent LU13 data. The blue open square indicates 
the central mass estimates for the cluster u Cen provided by Haggard et al. (2013), while the cyan open circle indicates the central mass 
hosted in the cluster NGC6388 as suggested by Lanzoni et al. (2013). All the models in this case have solar metallicities. 


To perform the simulations, we used the HiGPUs code 
(Capuzzo-Dolcetta et al. 2013), a direct IV-body integra¬ 
tor which fully exploits the advantages arising from parallel 
computing. HiGPUs allows very fast integration keeping a 
very high level of precision but, on the other hand, it does 
not implement any treatment of binary formation and close 
encounters and hence, it does not allow to follow the long¬ 
term evolution of the system. 

Moreover, the current version of HiGPUs does not con¬ 
tain any recipe for stellar evolution. To improve the compar¬ 
ison with our semi-statistical results, we developed a modi¬ 
fied version of HiGPUs in which we included the SSE pack¬ 
age. In the following, we will refer to this modified version 


as HiGPUsSE. We labelled in the following with SEn the 
simulations performed with HiGPUs, whereas we used SEy 
to label those carried out with HiGPUsSE. 

We set as gravitational softening t = 0.05 pc, thus al¬ 
lowing a reliable description of the evolution of the system as 
long as e is sufficiently smaller than the mean inter-particle 
distance, A, defined as (Gilbert 1968; Boily et al. 1999; Nel¬ 
son & Tremaine 1999): 

< 19) 

where F ~ 0.9 and n is the numerical density given by the 
stars which move in the innermost region of the cluster. 
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Figure 8. As in Figure 7, but here all the models have a value of the metallicity Z = 0.0004. 


A good estimate of A is provided by the lagrangian ra¬ 
dius which contains the 0.1% of the total GC mass. Hence, in 
the following we will use this typical radius to check whether 
our results provide a proper description of the dynamical 
evolution of the system. 

Since our methodology is not well suited to describe 
the long-term evolution of the cluster nucleus, we integrated 
these models up to ~ 10 2 crossing times, a time-scale suffi¬ 
ciently shorter than the relaxation time, but long enough to 
highlight and characterise the formation of an MSS. 

Stars’ masses are sampled according to a Salpeter or a 
Kroupa IMF. Moreover, to model the clusters we used either 
a Dehnen density profile or an uniform density profile, as 
we have done for the statistical work discussed above. As 


shown in Table 6, each model is labelled with the number 
of particles used to represent it, the kind of code used to 
evaluate dynamics and stellar evolution, the IMF used to 
sample stellar masses and the density distribution used to 
sample stars’ positions. 


4-0.1 Dynamical friction 

Direct JV-body simulations are powerful tools to follow the 
mass segregation of stars in stellar clusters, since in this case 
df comes out naturally as a consequence of the two-body in¬ 
teractions. On the other hand, the presence of a softening 
length that smooths the gravitational interaction may lead 
to less efficient energy exchanges among particles, altering in 
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Table 6. 


Parameters of the 7V-body simulations. 


Model ID 

code 

IMF 

p(r) 

7 

r GC 

(pc) 

Rgc 

(pc) 

Moo 

(10 3 M 0 ) 

Mac 
(10 3 M 0 ) 

N 

8 k_SEn_U_S 

HiGPUs 

S 

U 

0 

5 

5 

0.3 

0.3 

8024 

8 k_SEy_U_S 

HiGPUsSE 

S 

U 

0 

5 

5 

0.3 

0.3 

8024 

16k_SEn_U_S 

HiGPUs 

S 

u 

0 

5 

5 

5.1 

5.1 

16384 

16k_SEy_U_S 

HiGPUsSE 

S 

u 

0 

5 

5 

5.1 

5.1 

16384 

32k_SEn_U_S 

HiGPUs 

S 

u 

0 

5 

5 

11 

11 

32768 

32k_SEy_U_S 

HiGPUsSE 

S 

u 

0 

5 

5 

11 

11 

32768 

262k_SEn_D_S 

HiGPUs 

S 

D 

0 

0.7 

10 

100 

90 

262144 

262k_SEy_D_S 

HiGPUsSE 

S 

D 

0 

0.7 

10 

100 

90 

262144 

524k_SEn_D_K 

HiGPUs 

K 

D 

0 

1.6 

10 

500 

335 

524288 

524k SEy D K 

HiGPUsSE 

K 

D 

0 

1.6 

10 

500 

335 

524288 


Column 1: model name. Column 2: AT-body code used. Column 3: IMF used to sample the mass of each star: Salpeter (S) or Kroupa 
(K) mass function. Column 4: density profile used to sample stars’ positions: uniform (U) or Dehnen (D) density distribution. Column 
5: slope of the density profile. Column 6 : scale radius of the model in pc. Column 7: total radius of the cluster. Column 8 : mass of the 
cluster for r cu t —>■ oo. Column 9: mass of the cluster actually sampled. Column 10: number of particles used to represent the cluster. 


Table 5. 

Slope of the scaling relation between MSS masses and the host 
cluster masses. 


Model name 

a 

e a 

b 

£(> 

A1 

1.000 

0.002 

-2.23 

0.01 

A2 

0.996 

0.002 

-2.25 

0.01 

A3 

0.999 

0.001 

-2.045 

0.008 

A4 

0.998 

0.001 

-2.075 

0.006 

A5 

1.01 

0.01 

-4.233 

0.082 

A 6 

0.98 

0.04 

-4.86 

0.27 

B1 

1.000 

0.001 

- 2.202 

0.007 

B2 

0.999 

0.001 

-2.238 

0.009 

B3 

0.999 

0.001 

-2.024 

0.005 

B4 

1.001 

0.001 

-2.066 

0.007 

B5 

0.99 

0.01 

-4.160 

0.087 

B 6 

0.97 

0.03 

-4.77 

0.20 


Column 1: model name. Column 2-3: slope of the fitting function 
and relative error. Column 4-5: offset of the fitting function and 
relative error. 

some cases the orbital decay process. Due to this reason, we 
ran several test simulations of the same model at decreasing 
values of e. We find that the orbital evolution of the most 
massive stars is substantially unaltered when e < 0.05 pc. 
Hence, to obtain a good compromise between the computa¬ 
tional cost and the accuracy of the simulations, we set this 
limiting value as softening length. 

An accurate ./V-body modelling of the cluster would al¬ 
low to get the dynamical friction time-scale of stars. Hence, 
we tested the validity of our treatment for df comparing the 
theoretical estimates of tdf with the values obtained directly 
from the simulations. In particular, we selected in simulation 
262k_SEy_D_S those stars for which the orbital decay occurs 
over a time sufficiently greater than the time resolution of 
the simulation itself. This requirement is crucial in order to 
extrapolate the decay time of a stars in the IV-body model. 
Having detailed informations about the orbital evolution of 
these stars, we evaluated their df times, tdf.NB, comparing 
them with the theoretical estimates obtained through Equa¬ 
tion 7. In particular, Figure 9 shows the relative difference, 


£df, between the TV-body and the theoretical time-scale pro¬ 
vided by Equation 7: 

£df = 1 — tdf/tdfNB- (20) 

The maximum difference Edf between semi-analytic 
treatment and N-body simulations is ~ 20%. Extrapolat¬ 
ing the orbital parameters of stars is the main source of 
uncertainty for tdf.NB, due to the fact that our simulations 
have a poor time-resolution. 

The good agreement between the semi-analytical esti¬ 
mates and the simulations suggest that Equation 7 can be 
used to reproduce reliable sets of initial conditions for the 
cluster immediately after the completion of the mass seg¬ 
regation process, thus allowing to simulate only the subse¬ 
quent evolution of the cluster. 

It should be noted that a detailed simulation of the 
mass segregation process for a reliable cluster model with 
a mass above 10 5 M@ may require more than a month to 
be performed, even using the most advanced direct A-body 
codes and hardware facilities currently available. Hence, our 
approach can allow to gain a significant amount of time for 
simulating the GC evolution. 

As the mass segregation proceeds, heavy stars concen¬ 
trate toward the innermost region of the host cluster, lead¬ 
ing to a progressive increase of the density, which in turn 
enhances the probability to have close encounters among 
stars. 

In HiGPUs and HiGPUsSE, gravitational encounters 
are smoothed through the gravitational softening e. Hence, 
we stopped our simulations when the lagrangian radius 
which contains 0.1% of the total mass of the cluster became 
comparable to e, in order to reproduce in the most reliable 
way the early evolution of the host nucleus. 

For instance, looking at the time evolution of lagrangian 
radii shown in Figure 10 makes evident that the innermost 
lagrangian radius reaches a size comparable to the softening 
in f ~ 100 Myr for simulation 262k_SEn_D_S and t > 200 
Myr for simulation 524k_SEn_D_K. 
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Figure 9. Relative error, erjf = 1 — tdf/tdfNBi as a function 
of the star masses, between the df time-scale evaluated through 
Equation 1 and through the 7V-body simulation in the case of 
model 262k_SEy_D_S. 




Figure 10. Evolution of the lagrangian radii as a function of 
the time. Top panel refers to simulation 262k_SEn_D_S, whereas 
bottom panel refers to simulation 524k_SEn_D_K. In each panel, 
from bottom to top, each line represents the lagrangian radius 
containing 0.1 — 1 — 5 — 10 — 25 — 50 — 75% of the total mass of 
the cluster. 


4-0.2 Stellar evolution 

In this section we try to determine if stellar evolution can 
affect significantly the formation of an MSS. In order to do 
this, we will compare the simulations performed with HiG- 
PUs with those obtained with its modified version, HiG- 
PUsSE. On the other hand, it should be stressed that at 
the time, nor HiGPUs neither HiGPUsSE include any treat¬ 
ment for close interactions and do not account for NSs’ natal 
kicks. Hence, we limit our comparison to the time interval 
over which the MSS forms and without considering dynam¬ 
ical kicks suffered by NSs at their birth. 

Figure 11 shows the mass enclosed within mss for all 
the simulations performed, Mmss- In particular, we com¬ 
pared in each panel simulations performed with HiGPUs 
(labelled with SEn), in which stellar evolution is evaluated 
“a posteriori” with respect to the dynamical evolution of the 
cluster, and those performed through HiGPUsSE (labelled 
with SEy), in which stellar evolution is accounted along with 
the dynamical evolution of the system. 

Our estimate of the length scale of the host cluster is 
affected by a ~ 20% uncertainty, leading to a corresponding 
uncertainty on the value of mss (see Equations 9 and 10). 

As expected, the mass of MSSs in simulations labelled 
with SEn is slightly higher than in simulations denoted with 
SEy, since in the latter case mass lost by each star decreases 
the df efficiency, thus increasing its orbital decay time-scale. 
On the other hand, it is quite evident that the mass es¬ 
timates provided using both HiGPUs and HiGPUsSE are 
comparable within the error. Hence, stellar evolution seems 
to not affect significantly the early phase of formation of 
MSSs. 

Moreover, it should be noted that the scaling relations 
provided in Section 3.2 allow to evaluate MSS masses that 
differ few percent from those evaluated through the A-body 
simulations. 

In the next section, we will focus the attention on the 
formation of an MSS in models composed of 262k and 524k 
particles simulated through HiGPUs, which allows a much 
faster integration with respect to HiGPUsSE. 


4-0.3 Formation of an MSS in globular clusters 

In configuration 262k_SEn_D_S the orbital decay of the heavi¬ 
est stars leads to the formation of a dense sub-system, MSS, 
with a half-mass radius of r h = 0.06 pc, which extends up 
to 0.1 pc. 

As shown in Figure 13, the fraction of stars with ini¬ 
tial masses below 30 Mg that move within mss rapidly 
decreases in time both in simulations 262k_SEn_D_S and 
524k_SEn_D_K, leading to the formation of an MSS mainly 
composed of heavy stars. 

In particular, model 262k_SEn_D_S contains a population 
of 68 stars with masses above 30 Mq, for a total mass of 3163 
Mq, that drives the formation of an MSS over a time-scale 
of 20 Myr. Figure 12 shows the distribution of stellar types 
within the MSS in this case assuming an age of 2 Gyr, in or¬ 
der to compare this results with the semi-analytical results 
discussed in the previous section. The MSS composition is 
really similar to the theoretical results shown in Figure 4, 
with a dominant fraction of NSs and BHs and a small per¬ 
centage of MS stars and WDs. 
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Figure 11. Mass enclosed within tmss as a function of the time for all the simulations performed. Filled red circles represent simulations 
performed with HiGPUs, whereas open black circles represent simulations performed with HiGPUsSE. 


During the MSS assembly, heavy stars transfer energy 
to lighter particles, which move outward, leading the total 
energy of the MSS to decrease reaching negative values over 
a time £ ~ 65 Myr. During this process, the MSS contracts 
continuously, the most massive stars form the innermost core 
of the MSS while lighter stars move outward. Due to this, the 
long-term evolution of the MSS depend only on the gravita¬ 
tional encounters among the most massive stars, which may 
halt the MSS contraction through the formation of binary 
systems and close encounters. 

On the other hand, it has been widely shown that close 
encounters and binaries should not affect significantly the 
global properties of the MSS (Heggie & Giersz 2014; Fregeau 
et al. 2006). In particular, Morscher et al. (2015) showed that 


strong interactions among binaries and single stars should 
not alter significantly the fraction of stellar BHs. 

The energy provided by the binary hardening leads the 
cluster to expand by a factor of ~ 2 — 2.5, depending on 
the metallicity of the host (Trani et al. 2014). Using Equa¬ 
tions 9 and 10, and accounting for this expansion the values 
of tmss presented above still agrees with observational es¬ 
timates (van der Marel & Anderson 2010; Haggard et al. 
2013; Lanzoni et al. 2013). 

Assuming for this model a value of the metallicity Z = 
0.0004, we found from the IV-body simulation that the newly 
born MSS should have a mass Mmss = 500 ± 80 Mg. On 
the other hand, making use of the scaling relations provided 
above for a similar model, we found Mmss = 513± 15 Mg, a 
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Figure 12. As in Figure 4, but for simulation 262k_SEn_D_S 
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Figure 13. Fraction of stars with masses below 30 Mq which 
move within the size of the MSS, as a function of the time, in 
models 262k_SEn_D_S (solid line) and 524k_SEn_D_K (dotted line). 
In both the simulations the whole population of light stars is 
pulled out from the MSS. 

value compatible with the one obtained through the A-body 
simulation. 

The formation process of an MSS can be argued also 
looking at the evolution of the host cluster’s structure. For 
instance, we found that its surface density profile E(A) oc 
R~ a , with R < 0.1 pc the projected distance to the cluster 
centre, gets steeper as the MSS formation proceeds, and the 
slope passes from a = 1 to a ~ 1.30±0.05 by the end of the 
simulation. Moreover, we found an increasing slope in the 
projected velocity dispersion, <r(R) oc R^ 13 . Indeed, we found 
that /3 passes from a nearly 0 value to 0.12 ± 0.01 within 
R < 0.1 pc. This findings agree with observations, since 
the observed mass excesses are often associated to power- 
law surface density profiles or projected velocity dispersions 
within the inner region of the host clusters (van der Marel 
& Anderson 2010; Noyola et al. 2010; Kamann et al. 2014). 

Regarding configuration 524k_SEn_D_K, instead, the 
population of stars heavier than 30 Mq has an initial mass 
of M = 2.7 x 10 4 Mq and is composed of 561 particles. 

In this case, mass segregation process is slower than in 
the previous model, as expected by Equation 2 (see Fig¬ 
ure 10). Furthermore, it should be noted that the stellar 
mass loss time-scale in model 524k_SEn_D_K is smaller than 


the segregation time-scale. Due to this, stars will reach the 
GC centre after they have lost most of their initial mass 
and, therefore, the subsequent evolution of the MSS will be 
substantially dominated by binary formation and close en¬ 
counters. 

The concentration of the stars leads to an MSS with 
a half-mass radius 0.1 pc, slightly larger than in the case 
262k_SEn_D_S. It is worth noting that, even in this case, the 
MSS size estimates provided in Section 3.1 are comparable 
to the size achieved through the A-body representation. 

As in model 262k_SEn_D_S, the fraction of stars with 
initial masses below 30 Mq approaches a value close to 0 
within the simulated time, leaving an MSS composed only 
of the most massive stars. 

For this model, assuming a metallicity Z = 0.0004, we 
found that the MSS should have a mass Mmss — 1960 Mq, 
a value 7% smaller than the value estimated through the 
scaling relations provided above. 


5 DISCUSSION AND CONCLUSIONS 

In this paper we investigated how the global properties of 
star clusters can affect the formation process, in their cen¬ 
tres, of massive sub-systems (MSSs), using a semi-analytical 
method and a series of direct A-body simulations. 

Despite the methodology used here do not account for 
the effects related to the formation of binary and multiple 
systems, which may alter significantly the dynamics of stars 
in the cluster, it should noted that only hard binaries com¬ 
posed of compact objects (stellar BHs or NSs) may inject 
enough energy to halt the contraction (Chernoff & Huang 
1996; Heggie 1975; Fregeau et al. 2006). Moreover, recent 
works have shown that the formation of binaries should not 
deplete the population of compact stars within the cluster 
centre (Mackey et al. 2008; Heggie & Giersz 2014; Morscher 
et al. 2015). Hence, the formation of binaries within the 
cluster centre should not affect significantly our estimates 
of MSS masses and size, which represent reliable upper lim¬ 
its to the observed mass excesses, despite the simplicity of 
our approach. 

In the following, we summarise our main results: 

• in order to understand how the global properties of a 
star cluster affect the formation of an MSS, we investigated 
mass segregation in 168 models of star clusters with different 
initial density distribution, total mass, IMF and metallicity, 
aiming to quantify their role; 

• we found that the mass growth of an MSS is charac¬ 
terised by three distinct phases: i) a fast phase during which 
Mmss increases through orbital decay of heavy stars with 
small apocentres; ii) a second phase during which Mmss de¬ 
creases as a consequence of stellar mass loss; iii) a slow phase 
during which Mmss smoothly rises in corrispondence of the 
occasional orbital decay of stars that move in a peripheral 
region of the cluster; 

• though our results indicate that nor the spatial distri¬ 
bution neither the metallicity play a significant role in de¬ 
termining the final mass of an MSS, they are much more 
effective in determining which kind of stars populate it. In¬ 
deed, we found that MSSs formed in GCs characterised by 
an initial uniform density distribution should host a sig¬ 
nificant fraction of MS stars, while in the case of Dehnen 
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density profiles the MSSs are populated mainly by NSs and 
BHs. Furthermore, we found that low metal MSSs contains 
more BHs than NSs, unlike MSSs with solar metallicity; 

• we have shown that the population of stellar BHs tends 
to form a tight nucleus, which extends up to 0.1 pc, quite 
independently on the GC properties. On the other hand, the 
distribution of other stellar type is more influenced by the 
spatial distribution of the host; 

• we provided scaling relations connecting the masses of 
the MSSs and their host clusters, showing that they are 
in overall agreement with previous theoretical and observa¬ 
tional works. In particular, we have shown that these cor¬ 
relations allow to obtain a value for the central mass excess 
of clusters u> Cen in excellent agreement with the observed 
values provided by van der Marel & Anderson (2010); 

• in order to test our semi-analytical approach we per¬ 
formed a series of direct A r -body simulations using the code 
HiGPUs. Moreover, we used HiGPUsSE, a modified version 
of HiGPUs in which we included stellar evolution, in order 
to highlight its effects on the MSSs formation process; 

• comparing these simulations with our theoretical ap¬ 
proach, we demonstrated that the treatment for df used here 
is well suited to describe the orbital decay of massive stars. 
Moreover, comparing the results obtained with HiGPUs and 
those obtained through HiGPUsSE, we found that stellar 
evolution seems to not alter significantly the dynamical evo¬ 
lution of the cluster nucleus, where the MSS form; 

• we have shown that our treatment allows to predict 
quite precisely which kind of stars will populate the MSS, 
as shown through the comparison with the A r -body runs; 

• the global properties of the host cluster reflect the on¬ 
going formation of an MSS. For instance, we observed in 
262k_SEn_D_S a significant increase in the slope of the sur¬ 
face density profile of the cluster and in the slope of the 
velocity dispersion; 

• the agreement found between our semi-analytical esti¬ 
mates and the direct N-body simulations indicates that the 
procedure followed in this paper can be used to provide reli¬ 
able initial conditions for simulating the long-term evolution 
of GCs. This would allow to significantly shorten the compu¬ 
tational time needed to model the long-term evolution of a 
GC, since the completion of the mass segregation process in 
a typical GC would require computational times that exceed 
one month. 
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